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1 Introduction 

The theory of integration of the initial value problem 
for ordinary differential equations is well developed. 
It includes various and effective methods (j3] , [7j , for 
example) that are widely used [6]. A new approach 
for solving the Cauchy problem is introduced in this 
paper. It is based on approximation of the solution 
by a sequence of functions and deriving n nonlinear 
functional equations from the initial value problem. 
The procedure does not result in consecutive elimi- 
nation of low order terms of the local error when in- 
creasing the order of spectral approximation. Thus, 
presented discretization is different from known im- 
plicit Runge-Kutta methods. 

We employ the nonlinear approximation for con- 
structing a one-step explicit algorithm. Description 
of the algorithm is illustrated by an example on 
numerical integration of a large-scale nonlinear 
dynamical system. 

2 Discretization of the Problem 

'This is a corrected version of the paper "A spectral 
method for Cauchy problem solving" published in Proceedings 
of MCME International Conference, Problems in Modern Ap- 
plied Mathematics, N. Mastorakis, ed., WSES Press, 2000, 
pp. 227-232. Alterations were made in the explanatory part 
of the paper, whereas core mathematical constructions were 
copied exactly. Footnotes and an Appendix were added when 
correcting. 



Let vector Yj and vector functions Y and F have 
dimension N. We consider the problem 

Y'{T) = F{T, Y{T)), 

Y{Tj) =Yj, T G [Tj,Te] , Ti < Te, (1) 

supposing that the function F has properties that are 
necessary for the following constructions. 

Let n £ N, h = Te — Tj, a„ = h/Xnn, and A„„ is a 
parameter that will be defined below. Making use of 
the substitutions 

t = {T-Ti)/an, y{t)= Y{Ti + ant), (2) 

f{t,y) = anF{Tj + ant,Y{Tj + ant)) (3) 
we reduce the original problem to the following 

y'{t)=f{t,y{t)), y{0)=Yj, t g [0, A„„] . (4) 

We look for an approximation y{t) to prob- 

lem Q in the form 

n 

Vnit) = Yi + anot + Y^ {Pr^,t) (5) 

Snj{fin,t)= [ £nj{Pnt)dt (6) 

Jo 

that satisfies the initial conditions in Q. The un- 
known vector coefficients a„o, CLnj have dimension 
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N, j3n = 1, £nj{t) are the polynomials of exponents, 
which are orthogonal on the semi-axis with the weight 
function 1. The polynomials can be defined in the 
following recurrent way [1]: 



£nn{t) = exp(-nt), 



= (2n - l)exp( 
= (6„jexp(t) 



n 



l)t) — 2nexp(— nt), 



a 



nj 



(2i + l)(n + j)(n- j + 1), 



-nj 



(2j - l)2j(2i + 1), 
= Aj{n^ + f + n), 



dnj = (2j-l)(n-j)(n + j + l) 



j = n— l,n — 2,...,2. The recurrence relations also 
can be extended to j = 1, and this results in defi- 
nition of the polynomial £no{t) which has the zeros 
Ans) s = l,...,n, Xnn IS the maximum non-trivial 
zero. The polynomial Snoit) is orthogonal to the sys- 



tem of functions £n{t) = {£nj{t)}^=i, but it is not 
integrable on the semi-axis. 

Obviously, limt^oo^nit) = 0, and one can come 
to the conclusion that Sn{t) is complete in L2[0,oo). 
Hence, nice properties of the system on the semi-axis 
may serve a useful purpose in approximation of a 
continuous function on a finite interval if the system 
is completed by unity. Accordingly, the sequence of 
functions l,t, £n{t) may be employed for approxima- 
tion of a continuously differentiable function on an 
interval. Such a deduction leads to approximation 
([5]) , ([g]) , and the integral in the right-hand side of ^ 
is calculated by the formula 



£njiPnt)dt 




l=j+l 



nl 



Satisfying the initial condition for the derivative of 
the approximation in ^ and making use of 

the polynomial values £nj{0) = (— 1)""-^ one can im- 
mediately show that 



,{t) = Yi+fot + J2anj{Snji/3n,t) - (-1)' 



(7) 



where /o = /(0, Yj). 

To find the coefficients a 
discrete form 



Unions) — f{tnsi Unions)) 



in n7\j we consider the 



(8) 



of the ordinary differential equation in Q at the col- 
location points tns = ^ns/Pn, s = 1, ... ,72. We are 
going to reduce equations ([s]) to a system of func- 
tional equations, and we quote below a few more 
properties of polynomials Snjif) that are necessary 
for computations. 

One can observe that the system of functions 
{£n U £no){t) generates the Gauss-type quadratures 
for exponents on the semi-axis. This results in the 
discrete form of orthogonality of the polynomials : 



n 

E 

s=l 



Pns£nj{Xns)£nl{Xns) — + 0) 



V m=l / 



(9) 



j,l = 1,2, ... ,n, 5ji is the Kronecker delta, Xns and 
Pns are the abscissas and weights of the quadrature. 
Also, the formula 



Y^ Pns£nl{^ns) — 1/^ 



s=l 

takes place. It stands for the integral of £ni{t) on 
the semi-axis. It may be noticed that the formula 
for the weights in Q is similar to that one in j2] 
for the original Gauss quadrature. The Gauss rule 
for exponents also follows straightforward from the 
substitutions 



l-2exp(-A„ 



2p„sexp(-A„s), (10) 



were Zns and Wns are the abscissas and weights of the 
original Gauss quadrature on the interval [—1,1]. 

Substituting approximation ([T]) to equations ([s]), 
multiplying the equations by Pns£nl{^ns) , adding 
them, making use of the properties described, and in- 
verting the matrix in the left-hand side of the equality 



^V. S. Chelyshkov, a variant of spectral method 
in the theory of hydrodynamic stabiUty, Hydromechanics 
{Gtdromekhanikaj, N 68, 1994, pp. 105-109, (in Russian). 
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developed one can obtain the following equations 

n n 



where the coefficients in ( 15 ) are calculated as follows 



1=1 



s=l 



where 

J ns — fit {tns)), 

( (-1)'2, j^l, 

-1, j = I, I odd, 
3, j = I, I even. 

Equation (11) together with approximation Q at 
t = tns form the system of N x n functional equa- 
tions for the components of the vectors a„j. How- 
ever, we will not examine such a spectral form of the 
discretization, but return to the original variables. 
Substituting ( 11 ) to ([T]) we represent the approxi- 



mation in the form 

n 

y„(t)= Yj+f,Qno{t) + J2fnsQns{t), (12) 



s=l 



where the functions in the right hand side are 

Qnoit) = i-ir -2j2Snj{(3n,t) | , (13) 

n 

Qns{t) = 2pns l£nl ( Ans ) X 
1=1 

n 

Y^A,iSnj{Pn:t)-{-l)H \ . (14) 



Approximation (12) - (14) gives the opportunity to 



evaluate 2/„(t) for any current t if the values of 
Vnitris) are known. 



Equating (12) at the points t = tnp, p = 1, . . . ,n 
and making use of substitutions ([2]), ([3]), we finally 
obtain the discrete analogue of initial value problem 
([T]) in the form 

Tnp — Tj -\- Vnph^ Ynp — Y(Tnp), 

n 

Ynp — Yi + anpohF{Ti, y/) + /i^ ^npsF(Tnsi Yns")j 

s=l 

(15) 



tnp/ -^nnj 



'^npO — Qnoitnp) / ^nm '^nps — Qnsitnp) / 



(11) Although discretization (15) resembles an implicit 



Runge-Kutta method, the procedure developed is not 
subjected to the basic idea of the method. In fact, 
low order terms in the expansion of the approximate 
solution in the Taylor series in h are not eliminated 
when increasing n. 

If F{T, Y[T)) is a linear function of Y then the 
precise solution of problem ( 15 ) may be found |^ ; 
in the general case, it is a difficult problem that 
requires application of methods for solving nonlinear 
functional equations. We will now describe an 
explicit algorithm that provides approximation to 
the solution of the initial value problem. 



3 A Recurrence Algorithm 



It follows from the first formula in (10) that the ze- 
ros Xns of the polynomial <?no(i) cluster near t = 0, 
when n increases. This indicates that computation 
can originate in the lowest order of approximation 
and then can be extended from point to point by 
drawing in the polynomials of successively increasing 
degree. 

We introduce the recurrence index, k, k = 1, . . . ,n 
and substitute k for n in the previous construc- 
tions, where it is necessary. We consider succes- 
sively elongated subintervals [0, Xnk] and collocation 
points tks = tksin), s = 1, . . . , fc on them such that 
tkk = We meet the requirements of the dis- 

crete orthogonality for the points t^s on each of the 
subintervals in two steps. At first, we consider that 
f3k = /3fc(?i) for k < n and determine /3k by satisfying 
the condition Pktkk = Xkk, so (3k = Xkk/Kk- Next, 
we subject the choice of tks for s < fc to the condition 
I3ktks = Xks, and therefore tks = Xks/(3k- 

We suppose that the integration step of size h is 



sufficiently small and reduce approximation (12) to 
the explicit form 

k 

yk{t)= Yj+foQko{t) + J2gksQks{t), (16) 

s=l 



^If F{T, Y{T)) = F{T), Ti = 0, and ft = 1 then equality 
( 15 1 contains the quadrature rule that is exact for the inte- 



grands {l,exp(-jA„„T)}, j = 1, 



on the interval [0, 1]. 
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y'k-iitks), s = i,...,k-i, 

(17) 



Thus, we use extrapolation outside of the interval 
[0,tfc-i,fc-i] for s = k. 
A function, yQ{t), is required to initiate calculating 



sequence (16), (|T7|. The first two terms in approxi- 



mation ([7]) suggest themselves, and we define 

voit) =Yj + /o7^oo(^), 7^oo(^) = t. 

For A; = 1 we have 

Voitu) = 1/ + '7^oo(^ll)/o, 
9ii = f{tii,yo{tn)), 

y^{t) =Yi+ foTZioit) + QuT^uit), 

where 

nioit) = Qio{t), 7^n(^) = Qii(t). 
If n = 1, then y^(tii) also should be calculated. 



Making use of (13), ilAh and (16) we calculate the 



derivatives in (17), and for 1 < k <n the recurrence 



relations are as follows 

Iksr = T^'k-l,ri'tks), 

s = l,...,k — l, r = 0, . . . , k — 1; 

k-l 

y'k-liiks) = IksofQ + ^lksrOrr: S = 1, . . . , /c - 1, 



r=l 



fc-1 



Vk-liiks) — Yl + T^k-l,oitkk)f + T^k-l,ritkk)grry 

r=l 

9kk = f{tkk, Vk-li^kk)), 
k 

Vkit) =Yi+ foUkoit) + drrT^krit), 



r=l 



where 



TZkoit) = Qkoit) + gko{t) 
TlkAt) = Qkr{t), r = 1, . . . , /c - 1; 
Tlkkit) = Qkk{t), 



and 



fe-i 



Gkr{t) = 'YlksrQks{t), r = 0,...,k-l. 



If k = n, then yn{tnn) also should be calculated. 

Making use of substitutions ^ , ([s]) and of previous 
computations we finally obtain the discretization in 
the original variables 



K^o = hF{Ti, Yi), 



(18) 



p-i 



Kap = hF \Ti + Unph, Yi + ^InpsKns , (19) 



s=0 



YniTi + h)= Yi + Y^ 

^nns 



(20) 



s=0 



fJ'nps — T^p~l,si^ipp) / ^nni P — 1, . . . 



The algorithm developed requires n + 1 calculations 
of the right-hand side of problem ([T]). Formulas 
(18) - (20) are the Runge formulas, but the al- 



gorithm does not represent an explicit Runge-Kutta 
method. Rather, it can be characterized as the ex- 
plicit Euler method which residual is reduced by the 
spectral component of the approximation. 



Algorithm (18) - (20) may have favorable proper- 



ties of stability and monotonicity at the expense of 
precision. 



4 Examples 

First, we examine the simplest explicit algorithm. On 



putting n = 1 in ( 18 ) - (M we find 



Kw = hF{Ti, Yi), 
Kii = hF (Tj + h,Yi + Ki 



10) 

The algorithm has the same order as the Euler 
method, but requires two calculations of the right- 
hand side of the problem. The local error at T = 
Ti + h is 



Yi 



-0.0573 



dF ^dF\ ,2 ^,,3, 
dt dYjr 



T=Ti 



s=l 



(21) 

For the Euler method, the constant in (21) is equal 
to 0.5. For the differential equation y' = —7?/ with 
7 > 0, the spectral algorithm holds monotonicity if 
h < 1.7943/7, while for the Euler method h < I/7. 
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We elaborated codes for algorithm (18) - (20) 



Calculation in successively elongated subintervals in- 
volves exponential polynomials of increasing degree, 
and approximations of degree n—1 and n can be com- 
pared at the end of a step. The comparison was used 
for fulfillment certain conditions and, as the result, 
for implementation of adaptive step-size control. 

We examined some algorithms of degree n < 16, 
and the next example represents a test on numeri- 
cal integration of a large-scale nonlinear dynamical 
system. The system was extracted from an initial - 
boundary value problem for the Navier-Stokes equa- 
tions that describes evolution of two dimensional dis- 
turbances in the laminar boundary layer near a flat 
plate [1], [2]. Calculations were performed for the 
Reynolds number R = 10^. The dynamical system 
has = 1220 degrees of freedom, and the initial 
value problem was integrated over a long time in- 
terval that displays complete development of distur- 
bances. An example of simulation of two dimensional 
statistically steady flow for n = 16 is shown in Fig- 
ure 1 on a representative time interval. Similar 
graphs also were obtained by a Runge-Kutta method. 




Figure 1. Dependence of skin friction on time 

Since the explicit algorithm is the result of approxi- 
mation of the solution by a sequence of functions, the 
solution can be interpolated in the integration inter- 
val and extrapolated outside of it. This could lead to 
constructing algorithms with translation. The col- 
location points are distributed non-uniformly on the 
interval, and the maximum distance between them 
is proportional to A„„ — An-i,n- Thus, an algorithm 
may include the shift by /ia = (1 — A„_i_„/A„„)/i, in- 
terpolation and extrapolation of the solution at the 



shifted points, and calculation of the right hand side 
of the problem at the end point of partially trans- 
lated interval. 

5 Conclusion 

A spectral approach for integration of the initial 
value problem is developed in this paper. The 
approach results in an implicit procedure and 
explicit algorithm. The algorithm was tested, and 
it is comparable with classical explicit methods. In 
addition, the explicit algorithm may provide a trial 
value for solving stiff problems. 
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7 Appendix]^ 

Let A'^ = 1. We consider the stabihty function (j^ page 41), 
R„{z), of explicit algorithm (18 1 - (20 1 in the complex plane 
z. For n—1 and n = 2 we find 



Ri{z) ^ 1 + ^ + 0.5573052 , 

R2{z) « 1 + z + 0.5339542^ + 0.0988462^. 

s This indicates that convergence of R„{z) to is fairly slow, 
and the stability domain expands along the negative part of 
the axis ^(z) = 0. 

Let n = 1. For implicit procedure ( 15 1 we find 



Ri{z) = 



1 + Az 
1- (1-^)2' 



A = l/ln2 - 1 « 0.442695, the procedure is A-stable, and 
lim^^_^ Ri(z) ^ -0.794349. 



Alas, for n = 2 procedure ( 15 I is not A-stable. Thus, ap- 
proximation (|12[) - (|14[) primarily serves as the basis for con- 



structing the explicit algorithm, and discretization ( 15 1 cannot 
be recommended for solving stiff problems if n > 1. 

Below we outline two approaches that lead, respectively, to 
A- and L-stable procedures. 

First, we employ the same type of approximation, but sub- 
stitute special sequences of orthogonal exponential polynomials 
to it. 

Let polynomials of exponents £^k^\t) be orthogonal se- 
quences on the semi-axis with respect to the weight function 
(exp(-t))"(l - exp(-f))'' such that 



[n — k)\ da;"~* 



(a;"+"+'=(l-i)''+"-'=). 



The author is presently an employee of Lane College, TN, 
USA. 

*E. Hairer and G. Wanner, Solving Ordinary Differential 
Equations II, Springer, Berlin, 1991. 
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Also, let Ai"/' 



be the maximum non-trivial zero of the poly- 
nomial f^g''^''(t). 

Considering /3 = cja, we introduce parameter u} that sets 
the abscissa of the weight function maximum, t = ln(l -I- uj). 
Then, by subjecting choice of a„ to the requirements 



7„ < 1 and 7„ = max A 



nn 1 



(22) 



we introduce parameter 7„. In general, equality in (221 can be 
reached if an,/3n G K; inequality holds if a„,/3„ are selected 
from the set of whole numbers or from a set of rational numbers 
of interest. Accordingly, for 7,1 = 1 we define functions 

and, for 7„ < 1, 



Similarly, we define sequences and Z^^^'it). Choice 

of ujn sets distributions of the zeros of the functions under 
consideration, and Z^^\l) = 2^q"'(1) = 0. System of func- 
tions (^Z^ U-E,jq"''^ (t) is easier to construct, in particular 
for a„,/3„ e N. 

Solving the initial value problem, we follow the approach 
analogous to ([5|,([6|. Precisely, we approximate the right-hand 

side f{t, y{t)) = g{t) of the problem by the functions Z,^ (t) 
making use of calculation of g{t) on the interval [0, 1] at the 
zeros of Zl^"\t) and at t = 0. 

For Q„, /3„ £ N, iJn = 1, and n = 2 we obtain 

«2 = /32 = 6, 72 = ln((15-f \/l5)/7), 

git) ^ 5(0)e2o(t) + g{i^i)e2iit) + 3(1)622(4), 



= 1-^1/72, Ml =ln((8 + V15)/7), 
and the interpolating functions are 



620 (i) 



621 (t) 



15 + 



1 30 - 

1 e 

7 

15 + 22^15 



30 

y 



72* _|_ ^g-272t 



15-H5\/l5 



622 (t) = +\/T5 + 



15 



7 

22/15 



15 



7 

15\/l5 



-272* 



-272* 



7 7 

Completing construction of the procedure, we find the Butcher 
tableau 





1^1 
1 


^1 — <ih^ 
1 - g/72 


—rvi — s/72 
-r + {q + s)/72 


rux + {q + s)/72 
r - s/72 




1 - g/72 


-r + (g + s)/72 


r — s/72 



where g = (8 -f ^)/14, g = (8 - Vl5)/14, 
s= 3(1 Vl5)/4, s= 3(1- ^)/4. 
Further, we derive the stability function 



15, 



R2{z) 



1 + Az + Bz^ 
1- {I- A)z + Cz'^ 



J. C. Butcher, Numerical Methods for Ordinary 
ential Equations, John Wiley & Sons, Chichester, 2003. 



yl = l + (3-2ri.i)/(272'), 

B = si.i/(372) - 25/(2172) + 1/(27!), 

C = -qrui/-y2 + qr /■y2 + 1/(27!), 

limz_>_cx) R2{z) ~ 0.543836, and come to the conclusion that 
the procedure is A-stable. 



Next result represents reconstruction of approximation ( 15 1 
that leads to a different type of the stability function. 

The collocation point i; = is necessary to employ for initia- 
tion an explicit algorithm, but for an implicit discretization we 
can drop direct calculation of the derivative of the solution at 
the left end of the integration interval. Realizing this approach 
and making use of £nk{t) for constructing the approximation 
we obtain for n = 2 the tableau 



ln(2 + 



1^1 


qui + r/l32 


svx - r/P2 


1 


q + r/132 


s - r//32 




q + r/132 


s - r/132 



where j32 — In (3 4- -\/3), vi 
g = (3 - V3)/6, r = ^/3/6, s : 
We also find that 

R2{z) 



■ 1 - Mi//32, Ml 
g 4- 2r. 

1 + Az 



1- {1- A)z + Bz^' 
A = g/ii//32, B = r/ii//32^, and the procedure is L-stable. 



One may mention that discretization ( |15[ ) contains quadra- 
tures with weights of mixed signs, while the procedures outlined 
in this appendix hold the quadratures with positive weights. 
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